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Abstract 


As  part  of  the  Offshore  Energy  Environmental  Research  (OEER)/DRDC  Atlantic 
project  “Feasibility  of  a  Marine  Vibroseis  System  to  Minimize  Potential  Impacts  of 
Seismic  Surveying  on  Commercial  Marine  Invertebrates”,  DRDC  Atlantic  has  been 
investigating  the  possibility  and  technical  feasibility  of  using  coherent  broadband 
sources  with  complex  waveforms  (followed  by  matched-filtering)  to  replace  high  am¬ 
plitude  impulsive  sources  used  in  marine  seismic  surveys.  This  report  discusses  a 
MATLAB  implementation  of  an  acousto-elastic  wavenumber  algorithm  for  layered 
geoacoustic  media.  This  code  is  then  used  to  generate  broadband  time  series  for  dif¬ 
ferent  geological  configurations  and  incident  waveforms.  The  relative  performances  of 
different  impulsive  and  extended-time  broadband  sources  are  investigated.  As  well, 
the  effects  of  the  upper  air/water  interface  and  the  inclusion  of  attenuations  in  the 
elastic  layers  are  modelled. 


Resume 


Dans  le  cadre  d’un  projet  de  Fassociation  de  recherche  environnementale  en  energie 
extracotiere  (Offshore  Energy  Environmental  Research  -  OEER)  et  de  RDDC  At- 
lantique  portant  sur  la  faisabilite  d’un  systeme  vibrosismique  marin  visant  a  reduire 
au  minimum  les  repercussions  potentielles  des  etudes  sismiques  sur  les  invertebres 
marins  commerciaux,  RDDC  Atlantique  a  etudie  la  possibility  et  la  faisabilite  tech¬ 
nique  de  remplacer  les  sources  impulsives  a  amplitude  elevee  utilisees  dans  les  etudes 
seismiques  marines  par  les  sources  a  large  bande  coherentes  a  formes  d’ondes  com¬ 
plexes  (suivies  d’un  hltrage  adapte).  Le  present  rapport  porte  sur  F integration  d’un 
code  MATLAB  a  un  algorithme  de  nombre  d’ondes  acousto  elastique  destine  a  des 
milieux  geoacoustiques  en  couches.  Le  code  permettra  ensuite  de  generer  des  series 
chronologiques  a  large  bande  pour  differentes  configurations  geologiques  et  formes 
d’onde  incidentes.  On  a  etudie  les  rendements  relatifs  de  differentes  sources  a  large 
bande  impulsives  et  a  duree  etendue.  On  a  egalement  modelise  les  repercussions  de 
l’interface  superieure  entre  Fair  et  l’eau,  ainsi  que  l’integration  d’affaiblissements  aux 
couches  elastiques. 
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Executive  summary 


Numerical  modelling  of  coherent  broadband  pulses  for 
seismic  exploration 

John  Fawcett;  DRDC  Atlantic  TM  2011-236;  Defence  Research  and  Development 
Canada  -  Atlantic;  November  201 1 . 

Background:  Typically  in  marine  seismic  exploration,  large  amplitude  impulsive 
sources  such  as  airgun  arrays  are  used.  There  is  concern  as  to  the  effects  that  such 
sound  has  on  marine  life  including  crustaceans.  There  is  the  possibility  that  one 
could  used  lower  peak  amplitude  pulses  that  are  extended  over  time  and  frequency 
(for  example,  a  Chirp  pulse)  and  regain  the  signal-to-noise  ratio  through  coherent 
matched- filtering.  This  report  describes  an  implementation  of  an  acoustic/elastic 
model  and  examines  the  coherent  matched-filtering  of  Chirp  pulses  for  some  examples. 
As  well,  a  Ricker  wavelet  is  used  to  represent  a  more  impulsive  type  source. 

Principal  results:  A  MATLAB  code  was  written  to  allow  for  the  computation  of 
seismic  reflection  time  series  at  an  array  of  receivers.  The  reflection  time  series  for 
short-duration  or  extended-duration  incident  pulses  were  computed.  It  was  found  that 
an  extended-duration  pulse  (after  match-filtering)  could  yield  seismic  sections  similar 
to  that  of  a  short-duration  pulse  given  that  the  two  pulses  were  approximately  over  the 
same  frequency  band.  However,  the  polarity  information  which  might  be  observed  for 
the  reflection  of  an  impulsive  wavelet  is  lost  in  the  matched-filtered/envelope  display. 
For  higher-frequency  sources  the  echo  levels  of  deeper  reflectors  are  reduced  due  to 
attenuation.  Higher-bandwidth  sources  are  able  to  better  resolve  thinner  layers.  The 
direct  and  surface-reflected  pulses  can  significantly  complicate  the  seismic  section.  In 
the  case  of  match-filtering  the  sidelobes  from  these  strong  events  could  mask  weak 
echoes. 

Significance  of  results:  The  results  show  that  extended  (and  lower  peak  amplitude) 
pulses  may  be  used  in  seismic  exploration.  There  are,  of  course,  technical  issues  with 
constructing  a  coherent  source  with  ample  power  in  the  frequency  range  corresponding 
to  an  airgun  array.  In  the  case  that  the  coherent  source  is  in  a  higher  frequency  band  it 
is  expected  that  the  echoes  from  the  deeper  reflectors  will  suffer  greater  attenuation. 
Some  care  needs  to  be  taken  with  the  match-filtering  so  as  to  reduce  the  sidelobe 
levels.  It  is  hoped  that  the  MATLAB  code  developed  in  this  report  will  be  useful  for 
further  studies. 

Future  work:  In  future  work  we  would  like  to  continue  our  numerical  simulations 
of  seismic  scenarios.  In  addition,  the  collection  of  real  data  to  validate  the  numerical 
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Contexte  :  En  exploration  sismique  marine,  on  utilise  generalement  des  sources 
d’impulsion  a  grande  amplitude,  comme  des  reseaux  de  canons  a  air.  Toutefois,  il 
existe  certaines  preoccupations  quant  aux  effets  de  ceux-ci  sur  la  vie  marine,  y  com- 
pris  les  crustaces.  II  est  possible  d’utiliser  des  impulsions  a  amplitude  moins  elevee  qui 
sont  etendues  sur  le  plan  du  temps  et  des  frequences  (p.  ex.  une  impulsion  modulee 
en  frequence),  et  de  recuperer  le  rapport  “signal  sur  bruit”  au  moyen  d’un  filtrage 
adapte  coherent.  Dans  le  present  rapport,  on  decrit  la  rnise  en  oeuvre  d’un  modele 
acousto  elastique  et  etudie  le  filtrage  adapte  coherent  de  quelques  exemples  d’impul- 
sions  modulee  en  frequence.  En  outre,  on  utilise  l’ondelette  de  Ricker  pour  representer 
une  source  de  type  plus  impulsive. 

Resultats  :  On  a  developpe  un  code  MATLAB,  afin  de  permettre  de  calculer  les 
series  chronologiques  de  sismique  reflexion  d’un  reseau  de  recepteurs,  ainsi  que  les 
series  chronologiques  de  reflexion  d’impulsions  incidentes  de  courte  et  de  longue  duree. 
On  a  determine  qu’une  impulsion  de  longue  duree  (apres  un  hltrage  adapte)  peut 
produire  des  coupes  sismiques  similaires  a  cclles  d’une  impulsion  de  courte  duree,  car 
la  bande  de  frequences  des  deux  impulsions  est  semblable.  Toutefois,  les  donnees  de 
polarisation  observables  lors  de  la  reflexion  d’une  ondelette  impulsive  sont  perdues 
dans  l’affichage  d’enveloppe  ou  de  hltrage  adapte.  L’attenuation  reduit  les  niveaux 
d’echos  des  rehecteurs  plus  profonds  des  sources  a  haute  frequence.  Les  sources  a  plus 
grande  largeur  de  bande  peuvent  representer  plus  facilement  des  couches  minces.  Les 
impulsions  directes  et  rehechies  en  surface  peuvent  compliquer  considerablement  la 
coupe  sismique.  Dans  le  cas  d’un  hltrage  adapte,  les  lobes  lateraux  de  ces  fortes 
impulsions  pourraient  masquer  de  faibles  echos. 

Portee  :  Les  resultats  montrent  que  des  impulsions  etendues  (et  une  amplitude  plus 
basse)  peuvent  etre  utilisees  en  exploration  sismique.  Evidemment,  la  conception 
d’une  source  coherente  a  haute  puissance  dont  la  gamine  de  frequences  correspond  a 
celle  d’un  reseau  de  canons  a  air  pose  quelques  problcmes  techniques ;  si  la  gamme 
de  frequences  est  plus  elevee,  on  s’attend  a  ce  que  les  echos  des  rehecteurs  profonds 
subissent  une  attenuation  plus  forte.  II  faut  faire  preuve  d’une  certaine  rigueur  pen¬ 
dant  le  hltrage  adapte,  ahn  de  reduire  les  niveaux  de  lobes  lateraux.  On  espere  que 
le  code  MATLAB  decrit  dans  le  present  rapport  s’appliquera  a  d’autres  recherches. 
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Recherches  futures  :  On  souhaite  continuer  les  simulations  numeriques  de  scenarios 
sismiques  et  acquerir  des  donnees  reelles  pour  valider  le  rnodele  numerique. 
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1  INTRODUCTION 


Many  typical  marine  seismic  exploration  sources,  air-guns,  boomers,  etc  produce  a 
low  frequency  (e.g.  30  Hz  in  the  case  of  air  guns),  short-time  duration  pulse  which 
is  directed  into  the  seabed  below.  The  resulting  echo  time  series  yields  information 
about  the  structure  of  the  earth  below.  Instead  of  using  a  large-  peak  amplitude 
short  duration  pulse,  another  approach  ( [1]-  [3] )  is  to  use  an  extended  time-frequency 
pulse  such  as  a  Chirp  or  pseudo-noise  signals.  The  received  time  series  is  then  match- 
filtered  to  regain  the  time-resolution  and  signal-to-noise  ratio.  An  advantage  of  using 
an  extended  pulse,  is  that  its  peak  amplitude  can  be  lower  than  that  of  the  impulsive 
waveform  due  to  the  gains  achieved  by  the  matched-filtering. 

In  seismic  surveys  it  is  often  desired  to  obtain  echoes  from  reflectors  several  hundreds 
of  metres  deep  or  more.  One  of  the  limiting  factors  in  the  propagation  is  the  attenua¬ 
tion  within  the  layers.  This  attenuation  in  the  simplest  model  is  defined  as  a  dB  loss 
per  wavelength  so  that  the  effect  of  the  attenuation  increases  with  frequency  (see,  for 
example,  Ref.  4  for  a  discussion  of  various  attenuation  definitions).  Some  fundamen¬ 
tal  questions  which  can  be  addressed  through  modelling  are:  (1)  what  are  the  relative 
performances  of  using  a  low-frequency  impulsive  source  and  using  an  extended-time 
pulse  in  approximately  the  same  frequency  range  followed  by  matched  filtering  and 
(2)  what  is  the  degradation  in  performance  as  the  dominant  frequency  of  the  source 
increased  and  (3)  how  does  the  layer  resolution  depend  upon  the  frequency  content  of 
the  source.  In  addition,  modelling  can  be  used  to  study  signal  processing  issues  such 
as  minimizing  sidelobes  (from  the  matched  filtering)  and  beamforming  (although  we 
do  not  investigate  beamforming  in  this  report). 

The  numerical  model  that  is  developed  is  based  upon  a  wavenumber  integral  repre¬ 
sentation  of  wave  propagation [3].  It  is  formulated  for  a  single  frequency  but  then 
can  be  run  for  a  set  of  frequencies.  The  earth  model  consists  of  water  layer  (with 
or  without  an  upper  air  surface)  and  a  layered  elastic  bottom.  The  reflected  field  at 
each  horizontal  wavenumber  and  the  reflected  field  in  the  spatial  domain  is  computed 
by  performing  a  wavenumber  integral.  The  time  domain  response  is  computed  from 
Fourier  synthesis  of  the  frequency  responses.  Other  authors  [4-6]  have  developed  nu¬ 
merical  codes  for  this  type  of  problem.  The  code  used  in  this  paper  is  a  simple-to-use 
MATLAB  implementation  that  we  developed  on  the  basis  of  the  theory  outlined  in, 
for  example,  Refs. 6  and  7.  The  development  of  this  implementation  has  the  advan¬ 
tages  that  it  will  be  easy  to  adapt  in  the  future  for  different  modelling.  Due  to  the 
fact  that  the  code  is  written  in  MATLAB  also  allows  one  to  easily  include  MATLAB 
graphical  and  built-in  function  capabilities. 

We  first  describe  the  model  that  we  have  implemented  and  then  we  generate  sets  of 
time  series  as  a  function  of  receiver  horizontal  offset  for  different  short  and  extended 
pulse  types.  In  the  case  of  the  extended  pulse  types  (Chirp  pulses)  the  resulting  time 
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series  are  the  matched  filtered,  the  analytic  signal  computed  from  the  Hilbert  Trans¬ 
form,  and  the  amplitude  of  the  envelope  displayed.  Layered  geoacoustic  models  with 
and  without  attenuation  are  considered.  In  addition  our  model  is  formulated  to  easily 
compute  the  response  with  and  without  the  effects  of  the  upper  water/air  interface. 
In  addition,  the  direct  arrivals  can  also  be  excluded  so  as  to  better  emphasize  the 
reflection  echoes. 
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2  Acousto-elastic  modelling 


We  will  consider  a  cylindrical  coordinate  system  (r,  z,  6)  (z  decreases  for  increasing 
depth  into  the  Earth)  and  we  will  consider  a  compressional  potential  0  and  a  shear 
potential  i/j  [5-7]  where  the  particle  displacements  u  are  given  by 

u  =  V0 

u  =  V  x  V  x  (0,0,0).  (1) 

Within  a  constant-velocity  elastic  layer  the  functional  form  of  the  potentials  is  given 
by  (for  zero’th  azimuthal  order): 

0  =  J0(kxr)(A~exp(—i \J  ( oj/cp )2  —  k%z)  +  A+exp{i\J  {w/cp)2  —  k%z) 

0  =  J0(kxr)(B~exp(-i \J (u/cs)2  -  k%z)  +  B+exp[i\J (< x/cs )2  -  k%z)  (2) 

where  kx  is  a  particular  horizontal  wavenumber  and  Jo  is  the  Bessel  Function  of  zero’th 
order.  The  layer-potential  coefficients  A~,  A+ ,  B~ ,  and  B+  will  be  determined  by 
the  continuity  equations  imposed  ta  the  top  and  the  bottom  of  the  layer.  It  can  be 
seen  from  Eq.(2)  that  for  kx  >  lo/cp  or  kx  >  ui/cs  there  are  decaying  or  growing 
exponential  terms.  At  a  horizontal  interface  between  2  elastic  media,  the  boundary 
conditions  are  that 


< 

=  Ur 

ut 

=  U~ 

Trz 

Uz 

=  Tzz 

If,  for  example,  the  upper  layer  was  a  fluid,  then  there  would  be  no  shear  potential  in 
that  layer  and  the  first  boundary  condition  of  Eq.  (3)  would  be  dropped.  The  symbol 
r  denotes  the  stress  tensor  and  we  also  will  use  the  relations  that: 


7~rz 

Tzz 


,duz 


du , 

+  iA 

duz 


A V  0  +  2p 


dz 


(4) 

(5) 


Using  these  definitions  determines  the  4  boundary  conditions  in  terms  of  the  4  un¬ 
known  coefficients  above  and  the  4  coefficients  below  the  interface, 


Ai ) 

( Ah  \ 

A7 

—  B+ 

Aj+ 1 

Bt 

W+i 

BU 

\B7  / 

\  Bj+ 1  / 

(6) 
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The  expressions  for  the  elements  of  the  matrices  Rj  and  are  derived  in  the 
Appendix. 

The  coefficients  for  layer  j  as  defined  in  Eq.(6)  include  the  exponential  vertical  factors 
(since  we  evaluated  the  expressions  as  if  the  interface  was  at  *  =  0).  However,  a  layer 
has  a  vertical  thickness  Az.  The  values  of  the  coefficients  at  the  bottom  of  the  layer 
(oj)s  can  be  related  to  those  at  the  top  of  the  layer  ( dj)r  by 


(  exp(— iAz^p)  0.0  0.0 

0.0  exp(iAz7p)  0.0 

0.0  0.0  exp(— iA^7s) 

\  0.0  0.0  0.0 


0.0  \ 

0.0 

0.0 

exp(f A^7s)  / 


(%)t 


(7) 


The  expressions  and  7S  denote  the  vertical  square  root  factor  in  the  exponentials 
of  Eq.(2).  The  problem  with  Eq.(7)  arises  when  the  exponentials  grow  and  this  can 
be  the  case  for  A~  and  B~ ,  when  the  square  roots  in  the  exponentials  have  a  positive 
imaginary  part.  If  the  layer  coefficients  A+  and  B+  are  defined  in  terms  of  their 
values  at  the  top  of  the  layer  and  A~  and  B~  are  defined  in  terms  of  their  values  at 
the  bottom  of  the  layer,  then  the  resultant  equations  for  the  coefficients  will  contain 
only  decaying  exponentials. 


Let  us  consider  a  stack  of  N  elastic  layers  bounded  above  by  a  water  halfspace  and 
below  by  an  elastic  basement.  There  are  a  total  of  4N  +  1  (water)  +  2  (basement) 
unknowns.  In  the  water  we  assume  we  have  the  specified  downward  incident  field 
and  the  unknown  upward  compressional  coefficient.  In  the  basement  there  are  the 
unknown  downgoing  compressional  and  shear  coefficients.  For  N  layers  there  are 
N  elastic  interfaces  resulting  in  4N  equations.  At  the  water/elastic  interface  there 
are  3  equations,  resulting  in  a  total  of  4N+3  equations  in  4N+3  unknowns.  At  the 
water/elastic  seabed,  the  downward  compressional  coefficient  dp  is  specified.  Thus  we 
take  the  column  of  the  equation  matrix  in  the  water,  corresponding  to  the  downward 
wave,  multiply  it  by  —dp  and  consider  this  as  the  right-hand  side  in  a  system  of 
equations.  The  solution  of  the  system  of  equations  yields  a  vector  of  coefficients.  It 
is  the  coefficient  corresponding  to  the  upward  compressional  coefficient  in  the  water 
column  which  is  of  interest;  for  an  incident  coefficient  of  unity  this  yields  R(kx',uj ), 
the  reflection  coefficient  at  a  given  wavenumber  kx  and  at  a  frequency  oj.  Let  us 
suppose  the  incident  coefficient  for  a  point  source  is  given  by 


exp(i7 p(kx)zs) 

2*7p 

then  the  spatial  reflected  wavefield  in  the  water  column  is  given  by 


p(r,  zr) 


1 

27T 


Jo(kxr ) 


R(kx,  tn) 

2*7 P 


exp(i'yp(zs  +  zr))kxdkx. 


(8) 

(9) 
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This  integral  has  square  root  singularities  and  possibly  poles  along  the  real-axis.  In 
order  to  avoid  this  problem  we  displace  the  contour  of  integration  slightly  down  the 
imaginary  axis  in  the  complex-plane.  In  particular,  we  define 

kx{t )  =  t-  ie silent /tend)  (10) 


and 


7T  t 

dkx(t)  =  (1  —  ie - cos(7r - ))dt 


tend 


tend 


(ii) 


where  e  is  typically  small  and  tend  is  the  largest  value  of  t  considered.  We  typically  set 
tend  =  ot2nf/cw  where  a  is  a  factor  greater  than  one  and  cw  is  the  sound  speed  of  the 
water  column.  The  integral  of  Eq.(9)  is  numerically  discretized  using  the  trapezoidal 
rule  with  discrete  values  of  t,  {tj},j  =  1, ...,  Nt.  Considering  the  discretized  integrand 
of  Eq.(9)  without  the  Bessel  function  as  a  vector  v  (for  fixed  source  and  receiver 
depths)  and  considering  the  Bessel  function  has  a  Nr  x  Nt  matrix  B  where  Nr  are 
the  number  of  receiver  ranges,  then  the  integral  can  be  written  for  multiple  receiver 
ranges  as 

p(rj,zr,zs)  =  Bv.  (12) 


This  is  the  reflected  field  -  the  spectral  contribution  for  the  direct  field  exp(i^p\zr  — 
zs\)/(2i'yp)  can  added  into  v  or  the  direct  term 


exp(icu/c\/ r2  +  (zr  -  zs)2) / (An^/r2  +  (zr  -  zs )2)  (13) 


can  be  added  to  the  final  result. 


Thus  far  the  source  is  in  an  infinite  upper  water  halfspace.  In  reality,  this  halfspace 
is  bounded  above  by  the  water/air  interface.  To  include  this  effect,  the  effective  field 
incident  upon  the  bottom  now  becomes  in  the  wavenumber  representation, 

(exp(z7Pza)  -  exp(i7p(2/i  -  za)) 

®7p 

where  h  is  the  depth  of  the  water  column.  For  R(kx)  the  reflection  coefficient  in 
the  water  without  an  upper  interface,  the  new  effective  reflection  coefficient  with  the 
upper  interface  is  given  by 


R  —  R  exp(i2'yph)  R  +  Rexp{i2^ph)Rexp{i2^ph)R  +  ... 


R 

1  +  exp(i2^ph)R 


(15) 


The  modified  receiver  term  becomes 


exp(i7P£r)  -  exp(i7p(i7P(2/i  -  zr)))  (16) 

Thus  in  this  case  we  have  a  new  discrete  wavenumber  spectral  vector  v  but  the  large 
Bessel  function  matrix  B  of  Eq.(12)  is  the  same  as  for  the  no  upper  interface  case. 
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Layer 

cp(m/s ) 

cs(m/s) 

pig/cm*) 

ap[dB/ A) 

as(dB/ A) 

H(m) 

Water 

1500 

- 

1 

0.0 

0.0 

(100m) 

Layer  1 

1750 

400 

1.3 

0.0  (.3) 

0.0  (.6) 

50m 

Layer  2 

2100 

800 

1.8 

0.0  (.1) 

0.0  (.25) 

20m 

Layer  3 

2300 

1000 

2.0 

0.0  (.1) 

0.0  (.25) 

500m 

Basement 

3000 

1200 

2.4 

0.0  (.1) 

0.0  (.2) 

- 

Table  1:  The  geoacoustic  parameters  used  in  the  modelling  of  this  report 


The  final  outputs  from  this  modelling  program  are  two-dimensional  matrices  of  fre¬ 
quencies  and  receiver  ranges:  (a)  one  matrix  is  the  reflected  field  for  no  upper  interface 
and  (b)  one  matrix  includes  the  effects  of  the  upper  interface  and  the  direct  arrival. 
In  the  case  of  the  direct  (range  =  Rd)  and  surface-reflected  direct  arrival  (range  = 
Rs)  the  known  analytical  expression 

exp  (iuj/cwRD)  exp  (iu/cwRs),  .  . 

1  AnRD  4t tRs  ’  1  ’ 

is  used  with  cw  denoting  the  water  sound  speed.  These  matrices  can  then  be  processed 
by  other  programs  which  combine  them  with  specified  source  functions  and  produce 
range/time  series. 

3  Numerical  Examples 


For  the  numerical  examples  of  this  report,  the  geoacoustic  parameters  are  listed  above 
in  Table  1.  In  this  table  cp  denotes  the  compressional  sound  speed  of  the  layer,  cs,  the 
shear  speeds,  p  the  densities,  ap  and  as ,  the  compressional  and  shear  attenuations 
(if  used  in  the  computation)  and  finally  the  layer  thicknesses  H.  The  attenuations 
will  only  be  used  in  some  of  the  computations.  There  are  two  quite  thin  near-surface 
sediment  layers  followed  by  a  thick  layer,  terminating  in  a  reflector  at  570m  depth. 

For  the  first  computations  we  show  in  Figs.  1-3  the  resulting  transmission  loss  curves 
as  computed  by  the  method  of  this  report  and  by  OASES  [2]  for  frequencies  ranging 
from  0.25  to  32  Hz.  The  wavenumber  integrals  are  computed  out  to  the  maximum 
(10cu/1500,  0.2)  for  each  frequency  with  20000  discretization  points.  The  geoacous¬ 
tic  parameters  of  Table  I  with  the  attenuations  are  used.  For  these  computations 
e  =  0.001.  There  are  some  very  small  differences  between  the  2  curves  for  the  low 
frequencies  -  the  overall  agreement  of  all  curves  is  excellent  (it  should  be  noted  that 
the  agreement  is  so  good  that  the  curves  may  appear  indistinguishable  in  the  figures). 
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0.25  Hz 


Range(m) 


Figure  1:  Transmission  Loss  curves  OASES  (blue),  MATLAB  code  (red)  for  frequen¬ 
cies  of  0.25  to  1  Hz.  The  2  lines,  blue  and  red,  are  difficult  to  distinguish  due  to  their 
very  close  agreement. 


2.0  Hz 


Range(m) 


Figure  2:  Transmission  Loss  curves  OASES  (blue),  MATLAB  code  (red)  for  frequen¬ 
cies  of  2  to  8  Hz.  The  2  lines,  blue  and  red,  are  difficult  to  distinguish  due  to  their 
very  close  agreement. 
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12.0  Hz 


Range(m) 


Figure  3:  Transmission  Loss  curves  OASES  (blue),  MATLAB  code  (red)  for  frequen¬ 
cies  of  12  to  32  Hz.  The  2  lines,  blue  and  red,  are  difficult  to  distinguish  due  to  their 
very  close  agreement. 

We  now  compute  the  time-domain  responses  for  400  ranges  between  5  and  2000  m 
and  4096  frequencies  between  0  and  4095  x  0.125 Hz  (the  0  frequency  is  reset  to  0.05 
Hz).  The  source  and  receivers  are  90  m  and  80  m  above  the  seabed  (or  10m  and  20m 
depth  respectively  when  the  upper  wate/air  interface  is  included)  .  The  fine  frequency 
sampling  means  that  in  the  time  domain  we  can  compute  out  to  8  seconds  in  duration. 
A  value  of  e  =  0.004  was  used  to  offset  the  contour  in  these  computations.  In  Fig. 4 
we  show  the  direct  arrival  pulses  as  computed  at  the  horizontal  distance  of  505  m 
(trace  101)  for  the  three  different  pulse  types  considered  in  this  report,  the  first  pulse 
type  is  a  30-Hz  Ricker  pulse.  The  second  type  is  a  Chirp  pulse  with  length  T  of  2 
seconds  and  we  consider  first  a  frequency  sweep  between  10  and  70  Hz  and  then  a 
sweep  between  50  and  250  Hz.  These  Chirps  are  tapered  with  a  Tukey  window  (i.e. , 
split  cosine  bell)  of  fraction  0.5.  The  source  spectra  are  normalized  such  that  in  the 
time  domain  the  sum  of  squares  of  the  corresponding  signal  is  equal  to  unity. 

S(t)  =  sin(2nt(f0  +  ^r)),  0  <t<T  (18) 

It  can  be  seen  that  to  contain  the  same  amount  of  energy,  the  longer  Chirp  pulses 
have  a  peak  amplitude  of  less  than  20%  of  that  of  the  impulsive  Ricker  pulse.  The 
propagation  computations  are  such  that  in  the  frequency  domain  the  direct  arrival 
has  the  form  _cxp0^Au..ff)  jn  FjgS.  5.7  we  show  the  computed  seismic  sections  for 
the  three  pulse  types  (be  seismic  section,  we  simply  mean  the  time-horizontal  off  set 
image).  For  the  Chirp  pulses,  the  results  are  the  amplitude  of  the  envelope  of  the 
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2.5 


Time(sec) 


Figure  4:  Three  pulse  types:  (a)  Ricker  wavelet  (30Hz)  (h)  tapered  10-70  Hz  Chirp 
(c)  tapered  50-250  Hz  Chirp.  For  display  purposes,  they  have  all  been  multiplied  by 
a  factor  of  1.E6 

match- filtered  output.  For  the  matched  filtering  the  known  incident  pnlse  is  used  as 
the  replica.  These  sections  do  not  include  the  effects  of  the  upper  water/air  interface 
or  the  direct  arrival  which  would  further  complicate  the  time  history. 

The  resulting  seismic  sections  are  qualitatively  similar  to  each  other  and  have  sim¬ 
ilar  amplitudes,  although  the  broader  band  image  of  Fig.7  shows  some  more  detail. 
There  is  a  time  offset  with  the  matched-filtering  results  due  to  the  half-width  of  the 
template  signal  which  we  have  accounted  for.  Also,  for  the  higher-frequency/larger 
bandwidth  section  of  Fig.7,  there  is  an  arrival  at  the  larger  offsets  (starting  at  approx¬ 
imately  1000m)  which  is  not  observable  in  the  lower  frequency  results  of  Figs. 5  and  6. 
This  corresponds  to  the  combination  of  the  single  reflection  off  the  interface  at  50m 
depth  (bottom  of  the  first  sediment  layer)  and  the  multi-path  consisting  of  an  upper 
reflection  off  the  seabed/layer  1  interface  and  two  reflections  off  the  50m  interface. 
These  2  echoes  arrive  very  close  in  time  at  the  longer  ranges  and  at  lower  frequencies 
destructively  interfere  -  at  the  larger  bandwidth  they  start  to  become  resolvable. 

In  comparison  with  Fig.5,  we  show  in  Fig.  8  the  seismic  section  for  the  Ricker  wavelet 
when  the  direct  arrival  and  all  the  interactions  with  the  upper  water/air  interface  are 
modelled.  It  can  be  seen  that  the  interpretation  of  the  section  in  terms  of  arrivals 
from  reflectors  and  those  which  are  air-surface/seabed  multiples  is  difficult.  In  addi¬ 
tion,  the  relative  closeness  of  the  source  and  receiver  (10  and  20m  respectively)  to  the 
upper  surface  means  that  there  are  4  paths  associated  with  each  reflection,  i.e. ,  direct 
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Figure  5:  The  seismic  section  for  the  30-Hz  Ricker  wavelet  -  top  surface  and  direct 
arrival  not  included.  The  values  have  been  multiplied  by  1.E6 


0  500  1000  1500  2000 

Range(m) 


Figure  6:  The  seismic  section  for  the  Chirp  10-70  Hz  (2  seconds).  The  absolute  value 
of  the  envelope  of  the  match-filtered  output  is  shown  after  multiplication  by  1.E6. 
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0  500  1000  1500  2000 

Range(m) 


Figure  7:  The  seismic  section  for  the  Chirp  50-250  Hz  (2  seconds).  The  absolute 
value  of  the  envelope  of  the  match-filtered  output  is  shown  after  multiplication  by 
1.E6 

incidence/direct  receive,  one  surface  reflection  incidence/  direct  receive,  direct  inci¬ 
dence/one  surface  reflection  receive,  and  one  surface  reflection  incidence/one  surface 
reflection  receive. 

We  now  implement  the  attenuations  listed  in  Table  I.  The  resulting  reflection  seismic 
section  for  the  50-250  Hz  Chirp  is  shown  in  Fig. 9.  The  resulting  reflection  seismic 
section  is  similar  to  the  no  attenuation  sections  of  Fig.  7  but  with  many  arrivals  now 
significantly  lower  in  amplitude.  However,  the  arrival  from  the  reflector  at  570m 
depth  is  still  evident. 

Finally,  we  show  some  of  the  individual  time  series  (traces)  in  more  detail.  In  Fig. 10 
two  traces,  with  the  Ricker  and  the  2  Chirps  (after  matched  filtering)  are  shown  for 
a  near  hydrophone  (55  m  horizontal  separation)  and  the  farthest  hydrophone  (2000 
m  separation).  These  are  for  the  case  of  no  attenuation.  For  trace  11  of  Fig. 10,  the 
first  echo  is  from  the  seabed  and  this  is  evident  for  all  3  pulse  types.  The  second 
echo  corresponds  to  the  first  interface  at  50m  depth  combined  with  the  third  interface 
20m  further  in  depth.  This  interface  is  not  distinguishable  in  the  case  of  the  Ricker 
and  the  10-70  Hz  Chirp,  but  is  distinct  for  the  wider  bandwith  50-250  Hz  Chirp. 
The  echo  from  the  lowest  interface  at  570  m  depth  is  received  at  approximately  0.6 
seconds  and  is  evident  or  all  3  pulses.  At  the  far  hydrophone,  the  echo  from  the 
deepest  interface  arrives  first  and  the  seabed  reflection  is  at  about  1.35  seconds.  In 
Fig.  11  the  same  results  are  shown  but  now  with  attenuation  included  in  the  layers.  As 
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Figure  8:  The  seismic  section  for  the  Ricker  wavelet  with  direct  arrival  and  upper 
surface  interactions  included.  The  results  are  shown  after  multiplication  by  1.E6 


0  500  1000  1500  2000 


Range(m) 


Figure  9:  The  seismic  section  for  the  50-250  Hz  Chirp  (direct  arrival  and  upper  surface 
interactions  not  included).  The  attenuations  of  Table  1  are  considered.  The  results 
are  shown  after  multiplication  by  1.E6 
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Figure  10:  The  Ricker  wavelet  (blue)  and  the  two  matched- filtered  results  10-70 
Hz(green)  and  50-250  Hz(red)  for  the  hydrophones  at  55  and  2000  m  horizontal 
separation.  There  is  no  attenuation.  The  results  are  shown  after  multiplication  by 
1.E6 

would  be  expected  the  echoes  from  the  seabed  and  the  shallower  interfaces  are  largely 
unaffected  by  the  inclusion  of  attenuation.  For  the  50-250  Hz  pulse  the  reflection  from 
the  interface  at  570m  is  significantly  weaker. 

An  advantage  of  the  impulsive  Ricker  wavelet  is  that  the  phase  of  the  reflection 
is  visually  evident.  For  example,  for  Trace  11  the  polarity  of  the  reflected  Ricker 
pulses  from  the  seabed  and  the  layer  at  50m  depth  are  the  same  as  the  incident 
pulse,  whereas  the  reflection  of  Trace  400  at  about  1.4  seconds  has  approximately  the 
reverse  polarity.  In  the  case  of  the  match-filtered  results  the  phase  of  the  reflections 
are  not  observable. 


4  Discussion  of  Results 


We  have  written  a  MATLAB  program  for  the  computation  of  reflection  time  series 
in  a  horizontally  layered  geoacoustic  medium.  This  code  is  similar  in  its  formulation 
to  other  numerical  codes.  However,  because  it  is  written  in  MATLAB  it  is  easily 
incorporated  with  the  other  MATLAB  routines  and  graphics.  In  addition,  it  can  be 
easily  modified  in  the  future  to  meet  additional  modelling  requirements. 

Using  this  code,  the  time  series  for  an  array  of  receivers  and  for  different  pulses  were 
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Figure  11:  The  Ricker  wavelet  (blue)  and  the  two  matched- filtered  results  10-70 
Hz(green)  and  50-250  Hz(red)  for  the  hydrophones  at  55  and  2000  m  horizontal 
separation.  The  attenuations  have  been  included.  The  results  are  shown  after  multi¬ 
plication  by  1.E6 


computed.  In  the  case  of  long  extended  pulses,  the  outputs  were  match-filtered  with 
a  resulting  increase  in  output  amplitude  and  a  compression  in  time.  This  type  of 
modelling  allows  one  to  investigate  the  effects  of  centre,  frequency,  bandwidth,  and 
medium  attenuation.  In  addition,  one  can  also  consider  various  signal  processing 
issues  such  as  tapering  the  pulse  to  reduce  sidelobe  levels  in  the  matched-filtering. 
One  could  also  consider  more  advanced  deconvolution  techniques  to  improve  the 
resolution  of  the  resulting  seismic  sections. 
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Annex  A:  Details  of  the  scattering  matrix  R 


In  this  appendix  we  give  the  detailed  formulas  for  the  elements  of  the  scattering  matrix 
R  described  in  Section  2.  In  particular,  we  utilize  the  4  elastic  continuity  equations 
in  terms  of  downward  and  upward-going  compressional  and  shear  potentials.  Let 
use  take  the  upward  compressional  potential  (j)+  as  A+  J0(krr)  exp(i-^/(^)2  —  k2)-,  the 
downward  compressional  potential  as  A~J0(krr )  exp (— *\/(y)2  —  k2)  and  similarly  for 
the  shear  potentials  in  terms  of  B+  and  B~ .  Then  for  the  displacements  and  stresses 
we  obtain  (evaluated  for  2  =  0): 


d(j)  d2^ 

dr  drdz 

=  —  kr(A+  +  A~  +  irysB+  —  irysB~)J1(krr) 


(A.l) 

(A.2) 


The  kr.J\  (krr)  term,  arising  from  the  radial  derivative  of  the  potentials  is  common  to 
all  terms  for  ur  and  can  be  ignored.  For  the  vertical  displacement  uz  one  obtains 


dA  _  A  /A  ON 

dz  r  dr  dr 

=  (ilp(A+  -  A~)  +  k2(B+  +  B~))UKr)  (A. 4) 


For  the  stresses  we  have  that 


.  dur  duz 

=  R{—  + 


and 


CJ 7* 


dz  dr 


d2cj) 

=  M2  ^  „  + 


drdz 


dA  ,2  A) 

dz2dr  r  dr  ’ 


azz  =  A  V2</>  +  2  fi 


duz 

dz 


=  -\u2/c2p<t>  +  2A^  +  k2r^) 

=  ((-Ac o2/cl  -  2/j,A)(A+  +  A~)  +  2/j,k2i-fs(B^ 


(A.5) 


(A.6) 

=  —  krfi(2i'ypA+  —  2i'fpA+  +  (2k2  —  u2/c2)(B+  +  B~))Ji(krr)  (A. 7) 


(A.8) 

(A.9) 

B~))J0{krr)  (A. 10) 


From  these  expressions  we  can  construct  the  matrices  RX  and  Rj  discussed  in  Section 
2.  From  the  expressions  above,  we  obtain  for  R  (the  superscript  +  would  use  medium 
parameters  for  the  upper  medium  and  -,  medium  parameters  for  the  lower  medium) 


/  1 

*7 p 

2iq p 

\  -Xu2 /c2p  -  2^2p 


1 

~A p 

-2*7p 

—Xu2  / c2  -  2^1 


hs 

k2 

(2k2  -  u2/c2s) 
2jik2ils 


-As  \ 

k2 

(2k2  -  u2/c2s) 
-2nk2ils  / 


Here  we  have  used  the  same  notation  as  in  Section  2. 


(A.11) 
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